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In this paper, we demonstrate the efficiency of simulations via direct computation of the parti- 
tion function under various macroscopic conditions, such as different temperatures or volumes. The 
method can compute partition functions by flattening histograms, through the Wang-Landau recur- 
sive scheme, outside the energy space. This method offers a more general and flexible framework for 
handling various types of ensembles, especially the ones in which computation of the density of states 
is not convenient. It can be easily scaled to large systems, and it is flexible in incorporating Monte 
Carlo cluster algorithms or molecular dynamics. High efficiency is shown in simulating large Ising 
models, in finding ground states of simple protein models, and in studying the liquid-vapor phase 
transition of a simple fluid. The method is very simple to implement and we expect it to be efflcient 
in studying complex systems with rugged energy landscapes, e.g., biological macromolecules. 
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In recent years, methods for Monte Carlo (MC) sim- 
ulation have been dramatically improved over the tra- 
ditional Metropolis algorithm A large class of MC 
methods are those based on the flat energy histogram, 
such as the multicanonical ensemble method the en- 
tropic sampling method the density of states (DOS) 
method and the statistical temperature method 
In this study, we demonstrate the efficiency of an alterna- 
tive sampling method, which simultaneously and directly 
computes the partition function at various values of a cer- 
tain macroscopic variable, e.g., temperature T or volume 
V. Since one does not know the partition function in 
advance, the partition function at different values of a 
chosen variable is initially set to unity and continuously 
modified throughout the simulation until convergence. 

We first demonstrate the case of sampling based on a 
number of discrete values of temperature. In this case, a 
number of sampling temperatures are set over the tem- 
perature range of interest. Similar to the expanded en- 
semble method or the simulated tempering method Q, 
two types of MC moves are used: an energy move un- 
der a fixed temperature and a temperature move under a 
fixed energy. Before each MC step, a fixed probability is 
used to determine which type of move the system takes. 
For the energy move, the Metropolis algorithm is per- 
formed at the present (reciprocal) temperature (3. For 
the temperature move, another temperature is ran- 
domly chosen, and the following acceptance probability 
is used to accept the move: 



Acc(/3 ^ /?') = min <^ 1 



eM-l3E)/Z0 



(1) 



Here E is the present energy; Zp and Zpi are the val- 



ues of the estimated partition function at temperatures 
f3 and f3', respectively. The partition function is "esti- 
mated" because it is unknown in advance. After each 
MC step, the estimated partition function at the present 
temperature is multiplied by a factor / > 1 [3| ■ This can 
be written as, 



In Zp \nZf} + In /. 



(2) 
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Similar to the WL algorithm, it is shown that by re- 
peating the above procedure for a fixed /, the estimated 
partition function can eventually converge within certain 
fluctuations proportional to y'lnf 0, Q . Moreover, due 
to the frequently modified acceptance probability, the ad- 
ditional errors in the estimated partition function (due 
to violation of the detailed balance condition) are larger 
in a stage with a larger In/. Therefore, the value of 
In / should be gradually decreased to improve the accu- 
racy of the estimated partition function. In practice, the 
whole simulation is separated into several stages, each 
marked by a different value of In / [3| . In passing from 
one stage to the next. In / is modified to {In f)/n [4]. We 
use n = ^/TO in this study so that In / is decreased by 
an order of magnitude every two stages (the procedure 
for optimizing the In / of each intermediate stage will 
be given in a forthcoming paper [Sj]). At the end of the 
simulation, In / is reduced to a tiny number such that vi- 
olation of the detailed balance condition is negligible. For 
each / stage, if the simulation runs for sufficient number 
of steps, each temperature receives on average an equal 
number of visits, i.e., a flat temperature histogram is 
achieved. Here the term "temperature histogram" refers 
to the number of visits to each discrete temperature in- 
stead of to a temperature interval. The simulation is 
allowed to enter the next / stage when the histogram 
fluctuation falls below a cutoff percentage ■'2) . 

An alternative approach is to fix the number of simu- 
lation steps by C/^/ln / for an / stage. It can be shown 
that the two approaches are equivalent for sufficiently 
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long simulations [8|. The constant C can be estimated 
from a few initial / stages. The second approach ensures 
a better convergence for a stage with a smaller In/. 

In principle, any set of sampling temperatures of in- 
terest can be used. However, two consecutive tempera- 
tures must be close enough to allow sufficiently frequent 
temperature transitions. This requires a certain over- 
lap between the energy distributions of two neighbor- 
ing tem perature s. This condition can be expressed as 
AT - ^{AE^)/Cv T/VC^ , where Cy and y/{AE^ 
are the heat capacity and energy fluctuation at temper- 
ature T, respectively. Therefore, the number of sam- 
pling temperatures is roughly proportional to VTV (ex- 
cept around the critical region), where N is the system 
size. This feature is advantageous for larger systems, 
which is also a merit of the parallel tempering method , 
but the latter does not deliver the partition function 
quickly. 

The algorithm was first tested on the 256 x 256 
square lattice Ising model. A wide temperature range, 
T G [0,8], was simulated in a single simulation. Since 
the sampling temperature increment of an efficient 
simulation should be inversely related to the heat 
capacity as discussed above (nonuniform temperature 
setup is known to be advantageous for this 

large system, sampling temperatures were distributed 
based on the roughly estimated heat capacity (e.g., 
that from simulation of a smaller system). Accord- 
ingly, the entire temperature range was partitioned 
into 13 subranges. Sampling temperatures were lin- 
early distributed inside each subrange with a different 
increment. The temperature subranges and their incre- 
ments were (0.1, 1.0|0.1), (1.0, 1.8|0.04), (1.8, 2.0|0.02), 
(2.0,2.2|0.005), (2.2,2.25|0.0025), (2.25, 2.3|0.002), 
(2.3,2.35|0.005), (2.35, 2.5|0.01), (2.5, 2.7|0.02), 
(2.7,3.6|0.05), (3.6,5.0|0.07), (5.0, 6.0|0.1), and 
(6.0, 8.0|0.2). Here the notation for each subrange 
is (beginning temperature, ending temperature | incre- 
ment). In total, there were 218 sampling temperatures. 
Each time the probability of choosing temperature over 
energy moves was 0.1% (this number should be larger 
for smaller systems). The modification factor In/ was 
decreased from 1.0 to 10 the number of MC steps for 
stage / was 100/Vln/ sweeps, so the whole simulation 
took 7.2 X 10^ sweeps. Thermodynamic quantities at 
temperatures other than the sampled temperatures can 
be calculated using the multiple histogram method [ll|- 
Histograms from the last / stage were used. The exact 
results of the Ising model were also calculated using 
the method by Ferdinand and Fisher [l3|. The relative 
errors of the partition function, energy, entropy, and heat 
capacity were no larger than 0.00064%, 0.071%, 1.1%, 
and 3.9%, respectively. Fig. [1] shows the results for the 
partition function and heat capacity. For comparison, 
the WL algorithm was applied to the same system 
using 15 independent simulations, and the maximum 
relative errors of the free energy, energy, entropy, and 
heat capacity were 0.0008%, 0.09%, 1.2%, and 4.5%, 



respectively The simulation cost of the WL algo- 
rithm was 6.1 X 10^ sweeps However, the acceptance 
probabilities for energy moves can be precalculated to 
avoid expensive exponential computation in our case. 
The above simulation was finished in 10 hours on a 
single Intel Xeon processor (2.8 GHz). 
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FIG. 1: Results for the 256x256 Ising model. The upper panel 
shows the partition function as a function of temperature. 
The curve is shown for In Z per spin with the contribution 
of the two ground states subtracted. The lower panel shows 
the heat capacity per spin as a function of temperature. The 
relative errors are shown in the insets for both panels. 

Next, we introduce a variation of the above algorithm 
that tries to find the transition temperature automat- 
ically and to spend more effort sampling around that. 
This feature is desirable if the transition temperature is 
not roughly estimated in advance. This can be achieved 
by using a modified updating scheme, to let the system 
visit each temperature with a different frequency wp. In 
the acceptance probability Eq. ([T]), the values, Zp and 
Zp', of the estimated partition function are replaced by 
Zp/wfj and Zpi/wp', respectively, whereas the updating 
scheme Eq. ([2]) is changed to InZp ^ \n Zp + In f /wp. 
The temperature histogram is constructed in such a way 
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that the total number of visits to a particular tempera- 
ture P is now divided by its associated frequency Wf^ . To 
focus sampling around the transition temperature, the 
frequency W/3 can be specified as an increasing function 
of the heat capacity. Since the values of the heat capacity 
are unknown in advance, they are updated at the end of 
each / stage and are used in the next stage. The modi- 
fied algorithm was tested on the same 256 x 256 Ising sys- 
tem. The frequency W13 at temperature /3 was set as the 
square of the heat capacity per spin. Sampling temper- 
atures were uniformly distributed over the whole range, 
T e [0,8], with a fixed increment AT = 0.002. The 
probability of choosing temperature over energy moves 
was raised to 10%. The value of In / was lowered from 
1.0 to VW X 10~^ The simulation was kept running at 
each / stage until the temperature histogram fluctuation 
was lowered below 50%. The last stage was purposely 
extended to 5.0 x 10*^ MC sweeps to accumulate more 
statistical data. Totally, 9.8 x 10^ sweeps were used. The 
relative errors of the free energy, the energy, and the 
heat capacity were no larger than 0.000 45%, 0.055%, 
and 4.0%, respectively. 

It is also possible to realize rejection- free, hence 
more efficient, temperature transitions. First, the 
relative probability at each temperature (3i, Pi = 
exTp{—/3iE)/Z/s., is calculated for the present energy E. 
Next, the accumulated probability for each temperature, 
Qi = J2j<iPj/12j Pj' is also calculated, to form a se- 
ries of brackets, [Qi-i,Qi), i — 1,2,..., with Qq = 0. 
If a uniform random number r G [0, 1) falls in the ith 
bracket, (3i will be chosen as the next temperature. This 
type of temperature move is analogous to the heat bath 
algorithm for energy moves jl3i] . It is relatively expensive 
because of many exponential calculations. However, this 
expense is negligible if a more expensive non-Metropolis 
algorithm is used for the energy move. As an example, 
the Swendsen-Wang cluster algorithm [3] was used as 
the energy move on large two-dimensional Ising models. 
To improve the efficiency, the energy and temperature 
moves were merged in such a way that each energy move 
was immediately followed by a rejection-free temperature 
move. Simulations were performed on critical tempera- 
ture windows estimated by |T — Td ~ L~'''. Here = 1 is 
the critical exponent, and is the critical temperature. 
About 10—20 sampling temperatures were distributed in 
each window. Parameters and results are listed in Ta- 
ble m The efficiency is clear in terms of the number of 
simulation steps required to reach the desired accuracy. 

Molecular dynamics (MD) can be used as an energy 
move as well. In this case, the probability of tak- 
ing temperature over energy moves is 50%. Constant- 
temperature MD (a length-5 Nose-Hoover chain fl^ with 
force-scaling [l^) is used as a (potential-)energy move [H. 
The thermostat temperature Tq was set to be 0.5. The 
simulations were used to find ground states of AB pro- 
tein models Il7| . We were able to find all known ground 
states [5i.lT8l.[l9il20|. and several new ones with lower en- 
ergies. Table [H] lists the new ground-state energies, and 



TABLE I: Results for Lx L Ising models using the Swendsen- 
Wang cluster algorithm 14] as the energy move. Maximum 
relative errors were calculated by assuming the errors at the 
left boundary to be zeros. Here, T- and T+ define the tem- 
perature window, and AT defines the increment. 



L 


(r-,r+]AT) 


MC steps 


e(ln Z) 


e{Cv) 


64 


(2.0, 2.9 ] 0.1) 


0.7 X 10" 


4.0 X 10"" 


1.6% 


128 


(2.1,2.6 1 0.05) 


2.0 X 10" 


1.2 X 10"" 


1.1% 


256 


(2.2,2.42 ] 0.02) 


2.9 X 10" 


3.6 X 10"^ 


1.4% 


512 


(2.2,2.34 0.01) 


3.1 X 10® 


1.0 X 10"'' 


1.0% 


1024 


(2.24,2.30|0.005) 


3.1 X 10® 


6.9 X 10"** 


1.4% 



Fig- 121 shows the corresponding configurations. Compar- 
ing our results (for model I [13) with those from the 
statistical temperature method [5| , the new ground state 
of the two-dimensional (2D) 55mer, Fig. [^^a), has a dif- 
ferent topology in the two inner strands; the new ground 
state of the three-dimensional (3D) 55mer, Fig. [2I^c), has 
a more compact configuration. In both cases, our ground 
states have black-black clusters (strong attractions) that 
are more favorably packed with no exposed black beads. 




(c) (d) 

FIG. 2: Lowest-energy configurations of AB proteins (black, 
A; white, B). (a) 2D, 55mer, model I. (b) 3D, 55mer, model I. 
(c) 3D, 34mer, model II. (d) 3D, 55mer, model II. 
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TABLE II: Lowest energies of AB proteins with Fibonacci sequences. Results are compared with those from the annealing 
contour Monte Carlo (ACMC) , the energy landscape paving (ELP) the conformational space annealing (CSA) [20|], 
and the statistical temperature molecular dynamics (STMD) 5]. 



protein 


ACMC 


ELP 


CSA 


STMD 


This work 


2D, 55mer, model I 
3D, 55mer, model I 
3D, 34mer, model II 
3D, 55mer, model II 


-18.7407 

-94.0431 
-154.5050 


-42.438 
-92.746 
-172.696 


-18.9110 
-42.3418 
-97.7321 
-173.9803 


-18.9202 
-42.5789 


-19.2570 
-44.8765 
-98.3571 
-178.1339 



The WL-type algorithms have also been applied to 
Lennard-Jones simple liquid systems [2l| through com- 
puting the multidimensional DOS. Here, we demonstrate 
that the simulation can be carried out using volume, in- 
stead of temperature, as the sampling variable, where 
the temperature and particle number are held constant. 
Each volume move can be implemented as a change of 
the scale of the system. Therefore, it is convenient to 



adopt reduced coordinates s — r/vV. 
function is factorized to the ideal gas part ^^g, 



The partition 

Zig, and 

a potential part Zy, i.e., Z = ZigZy, where Zy = 
(l/F^)/dr^exp[-/3C/(r^)] = / ds* exp[-/3C/(s^; y)]. 
Thus, we can dynamically compute the potential part of 
the partition function Zy, instead of Z, in the acceptance 
probability Eq. ([T]). This method was used to study the 
liquid-vapor transition of a 108-particle Lennard-Jones 
system with half-box truncation and periodic boundary 
conditions. After the simulation, the Helmholtz free en- 
ergy can be obtained through F = Fjg — In Zy /f3, and the 
Gibbs free energy profile under pressure p can be derived 
through G = F +pV, at each sampling volume (or den- 
sity). For each simulation under a fixed temperature, the 
transition pressure was first determined by equalizing the 
two minima on the Gibbs free energy curve; the values 
of liquid density p+ and vapor density p_ were also de- 
termined correspondingly. Simulations were performed 
under different temperatures T G [0.85,1.20], with in- 
crement AT = 0.01. To accurately determine the po- 
sition of coexistence densities, the sampling density in- 
crements Ap were 0.002 and 0.0005 around the roughly 
estimated liquid and vapor coexistence densities, respec- 
tively, whereas the transition region was filled by a larger 
increment Ap = 0.005. Typically, about 300 volume sam- 
pling points were used in a single simulation. The com- 
puted vapor-liquid coexistence curve is shown in Fig. [31 
The relation p±- pc ^ a\Tc-T\±b\Tc-Tf (the critical 
exponent (3 — 0.3258 22}) was used to extrapolate the 
critical temperature Tc and the critical density pc based 
on the corresponding power-law regions. The estimated 
critical temperature Tc and critical density pc were 1.304 
and 0.315, respectively. The results for this small sys- 
tem are consistent with those of the infinite system (e.g., 
Tc = 1.3123 and pc = 0.3174 [H). 

In summary, we have demonstrated the efficiency of 
simulations via direct computation of the partition func- 
tion. The method has a range of advantages. An im- 
portant one is in the ground-state-oriented applications, 
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FIG. 3: Phase diagram for the 108-particle Lennard-Jones 
system. The empty circles are results of simulations, the solid 
line is from power-law fitting, and the solid circle represents 
the estimated critical point for this small system. 



such as in the protein folding problem, in which case 
the WL algorithm suffers from lack of efficient sampling 
around the ground state. This is because the location of 
the ground state, and hence the proper energy range over 
which the sampling should be performed, is not known 
in advance. The efficiency of the WL algorithm will be 
further reduced if the energy landscape in the last en ergy 
bin (near the ground state) is continuous and rugged [24| . 
By contrast, sampling in the temperature space does not 
require a priori information about the ground state and 
can sample the vicinity of the ground state with desired 
accuracy. 

Our method can be viewed as a generalization of the 
DOS-based WL algorithm [3| since the DOS is indeed the 
partition function of the microcanonical ensemble. In the 
case of canonical versus microcanonical ensembles, for ex- 
ample, the partition functions of them are related by an 
expression, Z{N,V,T) = g{N,V,E)exp{-(3E)dE, 
where Z{N, V, T) is the canonical partition function and 
g{N, V, E) is the density of states or microcanonical par- 
tition function. It is easy to see that, in the canonical 
ensemble, one can fix any pair of thermodynamic param- 
eters and change the third one for sampling, while in the 
microcanonical ensemble, it is hard to do so, e.g., one can- 
not fix N and E to change V . This indicates that there 
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are inherent advantages in performing simulations (such 
as flattening the histogram) outside the energy space. Wc 
thus expect the general framework to be more flexible in 
handling other types of ensembles, especially the ones in 



which computation of the DOS is not convenient. 
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